home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_4 / a4_6.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  2.9 KB  |  116 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 4.6 (Chebyshev Approximation).
  9. % Section    4.5,    Chebyshev Polynomials, Page 246
  10. echo on; clc; format short; hold off; clear
  11.  
  12. % Investigation of a Chebyshev polynomial approximations.
  13.  
  14. % n+1 are points needed to construct Pn(x).
  15.  
  16. % The abscissas are stored in  X.
  17.  
  18. % The ordinates are stored in  Y.
  19.  
  20. % The points are counted k=1,2,...,n+1.
  21.  
  22. % Remark. cheby.m and tch.m are used for Algorithm 4.6
  23.  
  24. pause % Press any key to construct the Chebyshev coefficients.
  25.  
  26. clc;
  27. % Chebyshev polynomial approximation Pn(x) of f(x) over [-1,1].
  28.  
  29. % Place the degree of approximation in  n.
  30.  
  31. % Place the left endpoint in  a.
  32.  
  33. % Place the right endpoint in b.
  34.  
  35. % Place the 'string expression' for f(x) in  fun.
  36.  
  37. n = 4;
  38. a = -1;
  39. b = 1;
  40.  
  41. fun = 'exp(x)';
  42.  
  43. [C,X,Y] = cheby(fun,n,a,b);
  44.  
  45. pause % Press any key to form the coefficient polynomials.
  46.  
  47. clc;
  48. format short;
  49. Mx1 = 'Construction of a Chebyshev approximation polynomial.';
  50. Mx2 = 'The abscissas are:';
  51. Mx3 = 'The ordinates are:';
  52. clc,echo off,diary output,...
  53. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(X),disp(Mx3),disp(Y),...
  54. diary off, echo on
  55.  
  56. pause % Press any key to graph the Chebyshev approximation.
  57.  
  58. clc; clg;
  59. h = (b-a)/150;
  60. X1 = a:h:b;
  61. x = X1;
  62. Y1 = eval(fun);
  63. P = tch(C,n,x,a,b);
  64. plot(X,Y,'or',X1,P,'-g',X1,Y1,'-r');...
  65. hold on;...
  66. plot([a b],[0 0],'b',[0 0],[-10000 10000],'b');...
  67. xlabel('x');...
  68. ylabel('y');...
  69. Mx1 = 'Comparison of exp(x) and P';...
  70. Mx2 = [Mx1,num2str(n),'(x)'];...
  71. title(Mx2);...
  72. grid;...
  73. hold off;...
  74. shg; pause    % Press any key to continue.
  75.  
  76. Mx1='The Chebyshev polynomial has been rearranged in ordinary polynomial form.';
  77. Mx2='This ordinary polynomial looks like:';
  78. Mx3='Pn(x) = c(1)x^n + c(2)x^(n-1) + ... + c(n)x + c(n+1)';
  79. Mx4 = 'The degree is  n = ';
  80. Mx5 = ', and the coefficients list  C  is:';
  81. clc,echo off, diary output,...
  82. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(''),...
  83. disp(Mx3),disp(''),disp([Mx4,num2str(n),Mx5]),disp(''),...
  84. for i=1:5:n+1, disp(C([i:min(i+4,n+1)])); end,...
  85. diary off, echo on
  86.  
  87. pause % Press any key to graph the error for the approximation.
  88.  
  89. clc; clg;
  90. Z = zeros(1,n+1);
  91. plot(X,Z,'or',X1,Y1-P,'-g');...
  92. hold on;...
  93. plot([a b],[0 0],'b',[0 0],[-10000 10000],'b');...
  94. xlabel('x');...
  95. ylabel('y');...
  96. Mx1 = ['The error: ',fun,' - P'];...
  97. Mx2 = [Mx1,num2str(n),'(x)'];...
  98. title(Mx2);...
  99. grid;...
  100. hold off;...
  101. shg; pause    % Press any key to continue.
  102.  
  103. clc; format short e;
  104. X = -1:0.2:1;
  105. x = X;
  106. Y = eval(fun);
  107. P = tch(C,n,x,a,b);
  108. points = [X;Y;P;Y-P];
  109. Mx1=['Chebyshev polynomial approximation of f(x) = ',fun];
  110. Mx2='     x(k)         f(x(k))      Pn(x(k))     error';
  111. clc,echo off,diary output,...
  112. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  113. diary off,echo on
  114.  
  115.  
  116.